Joint segmentation and registration of images for object detection

ABSTRACT

A method for joint segmentation and registration of a contour for detecting an object in image data includes initializing a segmentation, initializing a registration of the segmentation, updating the segmentation and the registration, outputting an updated segmentation and an updated registration, wherein the updated segmentation and updated registration delineates the object in the image data.

This application claims priority to U.S. Provisional Application Ser.No. 60/549,459, filed on Mar. 2, 2004, which is herein incorporated byreference in its entirety.

BACKGROUND OF THE INVENTION

1. Technical Field

The present invention relates to computer aided diagnosis applications,and more particularly to a system and method for joint segmentation andregistration of images for object detection, for example, lymph nodedetection.

2. Discussion of Related Art

Detection of specific anatomic structures in medical images or volumesis an important research problem in Computer Aided Diagnosis (CAD)applications. Detection methods include methods for segmentation foroutlining a target structure and registration in the presence ofmultiple images of the same structure or region. In medical imagingapplications multiple number of image volumes in which a structure ofinterest resides may be available. Different modality images of the sameregion may also be available in some applications. The challenge is thento make use of and relate the existing extra information from severalgiven image volumes. Methods for dealing with the detection problem inmultiple images have been proposed. In one class of methods, a targetregion is segmented separately in both images and the resultingboundaries or regions are registered. In another class of methods, aglobal or local registration of the two images is carried out and thetarget region is segmented using information coming from both images.

However, no known system or method exists for a combined region andboundary-based flow for coupled registration and segmentation.Therefore, a need exists for a system and method for joint segmentationand registration of images for object detection.

SUMMARY OF THE INVENTION

According to an embodiment of the present disclosure, a method for jointsegmentation and registration of a contour for detecting an object inimage data comprises initializing a segmentation, initializing aregistration of the segmentation, updating the segmentation and theregistration, outputting an updated segmentation and an updatedregistration, wherein the updated segmentation and updated registrationdelineates the object in the image data.

The segmentation is an ellipse and the registration is rigid. Updatingthe segmentation and registration further comprises accumulatinggradients for the segmentation evolution for a plurality of images inthe image data and a plurality of samples on the ellipse, updating theregistration for each of the plurality of images, and updating thesegmentation given the updated registration.

The method comprises determining parameters of each rigid registrationincluding a translation vector, a rotation matrix and a non-uniformscale matrix.

The segmentation is a surface of a three-dimensional structure and theregistration is non-rigid. The non-rigid registration is a deformationvector field over the surface of the three-dimensional object. Updatingthe segmentation and registration further comprises accumulatinggradients for surface evolution for each coordinate in the image volume,and accumulating gradients for the non-rigid registration evolution foreach coordinate in the image volume. The deformation vector field mapsthe surface in a first image onto the surface of in a second image. Thethree-dimensional surface is determined from an initialized surface seedof the segmentation in a first image slice of the image data and twoneighboring slices in the image data in which the surface seed isextended.

According to an embodiment of the present disclosure, a program storagedevice is provided readable by machine, tangibly embodying a program ofinstructions executable by the machine to perform method steps for jointsegmentation and registration of a contour for detecting an object inimage data. The method steps comprise initializing a segmentation,initializing a registration of the segmentation, updating thesegmentation and the registration, and outputting an updatedsegmentation and an updated registration, wherein the updatedsegmentation and updated registration delineates the object in the imagedata.

According to an embodiment of the present disclosure, a method for jointsegmentation and registration of a contour for detecting an object in animage volume comprises initializing a seed in a first slice of the imagevolume, initializing a signed distance function, converting the distancefunction into a three-dimensional surface seed, initializing a non-rigidregistration for the surface, updating the distance function and thenon-rigid registration, and outputting an updated surface and distancefunction mapping the surface in the first slice to the surface in asecond slice, the surface delineating the object in the image volume.

Updating the segmentation and registration further comprisesaccumulating gradients for surface evolution for each coordinate in theimage volume, and accumulating gradients for the non-rigid registrationevolution for each coordinate in the image volume.

The surface is determined from an initialized surface seed of thesegmentation in the first image slice of the image data and twoneighboring slices in the image data in which the surface seed isextended.

BRIEF DESCRIPTION OF THE DRAWINGS

Preferred embodiments of the present invention will be described belowin more detail, with reference to the accompanying drawings:

FIG. 1 is a flow chart of a method according to an embodiment of thepresent disclosure;

FIG. 2 is an illustration of a system according to an embodiment of thepresent disclosure;

FIGS. 3A and 3B are flow charts of a method for object detection usingelliptical contour with rigid registrations according to an embodimentof the present disclosure;

FIGS. 4A-4I is an illustration of a coupled registration of threesynthetic ellipses using elliptical contour with rigid registrationsaccording to an embodiment of the present disclosure;

FIG. 5 is a flow chart of a method for object detection using a generalcontour with non-rigid registration according to an embodiment of thepresent disclosure; and

FIG. 6 illustrates a transformation according to an embodiment of thepresent disclosure.

DETAILED DESCRIPTION OF PREFERRED EMBODIMENTS

A combined region and boundary-based system and method performs coupledregistration and segmentation for detecting an object an image volume,e.g., an MR or CT image volume. According to an embodiment of thepresent disclosure, a finite dimensional registration domain is used forobject detection. The finite dimensional registration domain evolvesparameters of an ellipse and its registration parameters for objectdetection. According to an embodiment of the present disclosure, coupledPDEs (Partial Differential Equations) are used to estimate a contour ina source image or a segmentation, and its non-rigid deformation onto atarget image. The non-rigid deformation corresponds to a registration ininfinite dimensions without using shape priors.

It is to be understood that the present invention may be implemented invarious forms of hardware, software, firmware, special purposeprocessors, or a combination thereof. In one embodiment, the presentinvention may be implemented in software as an application programtangibly embodied on a program storage device. The application programmay be uploaded to, and executed by, a machine comprising any suitablearchitecture.

Referring to FIG. 2, according to an embodiment of the presentinvention, a computer system 201 for implementing a method for jointsegmentation and registration can comprise, inter alia, a centralprocessing unit (CPU) 202, a memory 203 and an input/output (I/O)interface 204. The computer system 201 is generally coupled through theI/O interface 204 to a display 205 and various input devices 206 such asa mouse and keyboard. The support circuits can include circuits such ascache, power supplies, clock circuits, and a communications bus. Thememory 203 can include random access memory (RAM), read only memory(ROM), disk drive, tape drive, etc., or a combination thereof. Thepresent invention can be implemented as a routine 207 that is stored inmemory 203 and executed by the CPU 202 to process the signal from thesignal source 208. As such, the computer system 201 is a general-purposecomputer system that becomes a specific purpose computer system whenexecuting the routine 207 of the present invention.

The computer platform 201 also includes an operating system andmicroinstruction code. The various processes and functions describedherein may either be part of the microinstruction code or part of theapplication program (or a combination thereof), which is executed viathe operating system. In addition, various other peripheral devices maybe connected to the computer platform such as an additional data storagedevice and a printing device.

It is to be further understood that, because some of the constituentsystem components and method steps depicted in the accompanying figuresmay be implemented in software, the actual connections between thesystem components (or the process steps) may differ depending upon themanner in which the present invention is programmed. Given the teachingsof the present invention provided herein, one of ordinary skill in therelated art will be able to contemplate these and similarimplementations or configurations of the present invention.

Elliptical Contour with Rigid Registrations Referring to the finitedimensional registration domain is used for object detection; given afinite number of images I_(i):Ω_(i)→

², i=1, . . . m, e.g., m=2 for illustration purposes; assume a commoncontour C lives on an independent domain Ω, and let C₁=g₁(C), C₂=g₂(C).Then the transformations that deform the common contour C onto contoursin image domains Ω_(i) are denoted by g_(i):

²→

², each of which is an element of a finite dimensional group G withparameters w_(i), . . . , w_(n). The goal is to find a closed curve CεΩthat deforms on the common domain whereas a curve C_(i)εΩ_(i)corresponding to the mapping C_(i)=g_(i)(C) deforms on the ith imagedomain Ω_(i). All contour parameters propagate according to a genericregion-based energy functional defined over both image domains asfollows: $\begin{matrix}\begin{matrix}{{E\left( {C,g_{1},g_{2}} \right)} = {{\int_{C_{1\quad i\quad n}}^{\quad}{{f^{1}\left( x_{1} \right)}\quad{\mathbb{d}x_{1}}}} + {\int_{C_{2\quad i\quad n}}^{\quad}{{f^{2}\left( x_{2} \right)}\quad{\mathbb{d}x_{2}}}} + {\int_{C}\quad{\mathbb{d}s}}}} \\{= {{\int_{C_{\quad{i\quad n}}}^{\quad}{{f^{1}\left( {g_{1}(x)} \right)}{g_{1}^{\prime}}\quad{\mathbb{d}x}}} +}} \\{{\int_{C_{\quad{i\quad n}}}^{\quad}{{f^{2}\left( {g_{2}(x)} \right)}{g_{2}^{\prime}}\quad{\mathbb{d}x}}} + {\int_{C}\quad{\mathbb{d}s}}}\end{matrix} & (1)\end{matrix}$

-   -   where f₁=f_(in) ¹−f_(out) ¹, and f_(in) ¹ and f_(out) ¹ are the        region descriptors inside and outside the transformed contour        g₁(C) respectively (same comments apply for f_(in) ² and f_(out)        ² for C₂). For example, a piecewise constant model can be        utilized by choosing        f_(i)=(I_(i)−mean_(in))²−(I_(i)−mean_(out))². A regularization        on the unknown contour C is typically included as given by the        last term above.

The evolution of the contour C can be given by: $\begin{matrix}{\frac{\partial C}{\partial t} = {\left\lbrack {{{f^{1}\left( {g_{1}(x)} \right)}{g_{1}^{\prime}}} + {{f^{2}\left( {g_{2}(x)} \right)}{g_{2}^{\prime}}} - \kappa} \right\rbrack N}} & (2)\end{matrix}$

-   -   where N denotes the unit normal to C, and k is the curvature of        the contour. The evolution of the first registration g₁ and the        second registration g₂ can be given by: $\begin{matrix}        {\frac{\partial g_{1}}{\partial t} = {\int_{C}\quad{{f^{1}\left( {g_{1}(x)} \right)}\left\langle {\frac{\partial{g_{1}(x)}}{\partial w_{i}},{\left( \left( g_{1}^{\prime} \right)^{- 1} \right)^{T}{g_{1}^{\prime}}N}} \right\rangle{\mathbb{d}s}}}} & (3) \\        {\frac{\partial g_{2}}{\partial t} = {\int_{C}\quad{{f^{2}\left( {g_{2}(x)} \right)}\left\langle {\frac{\partial{g_{2}(x)}}{\partial w_{i}},{\left( \left( g_{2}^{\prime} \right)^{- 1} \right)^{T}{g_{2}^{\prime}}N}} \right\rangle{\mathbb{d}s}}}} & (4)        \end{matrix}$

Referring now to an ellipse-segmenting flow with region-based energy,for specific applications in medical image analysis, the structures ofinterest may be approximated by an elliptical contour ε.

Incorporation of such a priori information into the evolution equationwill result in reduced degrees of freedom of the evolving curve andhence lead to a more robust estimation of the contour position and itsregistration.

A point on a 2D ellipse parameterized by pε[0,2π) is given by${p\begin{pmatrix}{a\quad{\cos(p)}} \\{b\quad{\sin(p)}}\end{pmatrix}},$and using its translation vector ${t = \begin{pmatrix}d \\e\end{pmatrix}},$and rotation matrix ${R = \begin{pmatrix}{\cos\left( \theta_{e} \right)} & {\sin\left( \theta_{e} \right)} \\{- {\sin\left( \theta_{e} \right)}} & {\cos\left( \theta_{e} \right)}\end{pmatrix}},$the parameterization of the elliptical contour in 2D is: $\begin{matrix}{{\varepsilon(p)} = {{{a\begin{pmatrix}{\cos\left( \theta_{e} \right)} \\{- {\sin\left( \theta_{e} \right)}}\end{pmatrix}}{\cos(p)}} + {{b\begin{pmatrix}{\sin\left( \theta_{e} \right)} \\{\cos\left( \theta_{e} \right)}\end{pmatrix}}{\sin(p)}} + \begin{pmatrix}d \\e\end{pmatrix}}} & (5)\end{matrix}$

Utilizing this parameterization, the variation of the energy in Eq. (1)with respect to ellipse parameters λ^(i)ε{a,b,d,e,θ_(e)}, i=1, . . . , 5yields the gradient flows: $\begin{matrix}{\frac{\mathbb{d}\lambda^{i}}{\mathbb{d}t} = {{\int_{\varepsilon}^{\quad}{{f^{1}\left( {g_{1}(x)} \right)}\left\langle {\frac{\partial{g_{1}(x)}}{\partial\lambda^{i}},N} \right\rangle\quad{g_{1}^{\prime}}{\mathbb{d}p}}} + {\int_{\varepsilon}^{\quad}{{f^{2}\left( {g_{2}(x)} \right)}\left\langle {\frac{\partial{g_{2}(x)}}{\partial\lambda^{i}},N} \right\rangle\quad{g_{2}^{\prime}}{\mathbb{d}p}}}}} & (6)\end{matrix}$for an evolution of the ellipse with this energy. ∂_(ε)/∂λ^(i) are givenfor each parameter by: $\begin{matrix}{\frac{\partial{\varepsilon(p)}}{\partial a} = {\begin{pmatrix}{\cos\left( \theta_{e} \right)} \\{- {\sin\left( \theta_{e} \right)}}\end{pmatrix}{\cos(p)}}} & (7) \\{\frac{\partial{\varepsilon(p)}}{\partial b} = {{b\begin{pmatrix}{\sin\left( \theta_{e} \right)} \\{\cos\left( \theta_{e} \right)}\end{pmatrix}}{\sin(p)}}} & (8) \\{{\frac{\partial{\varepsilon(p)}}{\partial d} = \begin{pmatrix}1 \\0\end{pmatrix}},{\frac{\partial{\varepsilon(p)}}{\partial e} = \begin{pmatrix}0 \\1\end{pmatrix}}} & (9) \\{\frac{\partial{\varepsilon(p)}}{\partial\theta_{e}} = \begin{pmatrix}{{{- a}\quad{\sin\left( \theta_{e} \right)}{\cos(p)}} + {b\quad{\cos\left( \theta_{e} \right)}{\sin(p)}}} \\{{{- a}\quad{\cos\left( \theta_{e} \right)}{\cos(p)}} - {b\quad{\sin\left( \theta_{e} \right)}{\sin(p)}}}\end{pmatrix}} & (10)\end{matrix}$

In Eq. (6), ${N(p)} = \begin{pmatrix}{{a\quad{\sin\left( \theta_{e} \right)}{\sin(p)}} + {b\quad{\cos\left( \theta_{e} \right)}{\cos(p)}}} \\{{a\quad{\cos\left( \theta_{e} \right)}{\sin(p)}} - {b\quad{\sin\left( \theta_{e} \right)}{\cos(p)}}}\end{pmatrix}$is the normal vector of the ellipse.

For each of the rigid registrations g_(i), a translation vector T, arotation matrix R, and a non-uniform scale matrix S are used. Thus:${\frac{\partial g_{k}}{\partial T^{x}} = \begin{pmatrix}1 \\0\end{pmatrix}},{\frac{\partial g_{k}}{\partial T^{y}} = \begin{pmatrix}0 \\1\end{pmatrix}},{\frac{\partial{g_{k}(x)}}{\partial\theta} = \begin{pmatrix}{- {\sin\left( \theta_{e} \right)}} & {\cos\left( \theta_{e} \right)} \\{- {\cos\left( \theta_{e} \right)}} & {- {\sin\left( \theta_{e} \right)}}\end{pmatrix}}$${\frac{\partial{g_{k}(x)}}{\partial S_{x}} = \left( {\begin{matrix}1 \\0\end{matrix}\begin{matrix}0 \\0\end{matrix}} \right)},{\frac{\partial{g_{k}(x)}}{\partial S_{y}} = \left( {\begin{matrix}0 \\0\end{matrix}\begin{matrix}0 \\1\end{matrix}} \right)}$

The evolution of the rigid registrations of the ellipse can be carriedout by Eqs. (3-4), or it can also be incorporated into the parametricellipse model to re-derive the evolutions of the rotation, translationand scale parameters from one image into the other. Eqs. (3-4) may beused to evolve the ellipse and the registration parameters.

Referring to FIG. 3A, an ellipse is initialized 301, for example, bymanually setting a seed position with present initial parameters for a,b, and θ. Rigid registrations are initialized 302 to the identitytransformations, e.g., a translation vector T=0, a rotation matrixR=Identity, and a non-uniform scale matrix S=Identity. Ellipseparameters and registration parameters are updated 303 for a time periodof interest. The ellipse is updated with parameters λ^(i) andregistrations are updated for g_(i)(TotalTime), yielding the ellipse onimages i=1, . . . , m 304.

Referring to FIG. 3B, the following flow may be used to updating asegmentation and registration 303. Δλ^(i)(t) = 0,j = 1, . . . , n(number of parameters for ellipse) For i = 1: m (number of images) (305)Δg_(i)(t) = 0 For p = 0:Δp:nSamples (loop around the ellipse) (306) x =coordinate of sample p //Accumulate gradients, see right portion of Eq.6 (307)${{\Delta\lambda}^{i}(t)}+={{f^{i}\left( {g_{i}(x)} \right)}\left\langle {\frac{\partial{g_{i}(x)}}{\partial\lambda^{i}},{{{g^{\prime}}_{i}}N}} \right\rangle{\Delta p}}$//Accumulate gradients, see right portion of Eqs. 3-4 (307)${{\Delta g}_{i}(t)}+={{f^{i}\left( {g_{i}(x)} \right)}\left\langle {\frac{\partial{g_{i}(x)}}{\partial w^{k}},{\left( \left( {g^{\prime}}_{i} \right)^{- 1} \right)^{T}{{g^{\prime}}_{i}}N}} \right\rangle{\Delta p}}$End //update the registration parameteres for each image (308) g_(i)(t +Δt) = g_(i)(t) + ΔtΔg_(i)(t) //discretized Eqs. 3-4 End //update thesegmentation of the ellipse (309) λ^(i)(t + Δt) = λ^(i)(t) + ΔtΔλ^(i)(t)//discretized Eq. 6

For ellipsoid flows, an extension of the flows in 2D to 3D is possibleby modifying the parameterization from an ellipse to an ellipsoid asfollows: ${\varepsilon\left( {s,p} \right)} = {{R\begin{pmatrix}{a\quad\cos\quad s\quad\sin\quad p} \\{b\quad\sin\quad s\quad\sin\quad p} \\{c\quad\cos\quad p}\end{pmatrix}} + d}$where sε[0,2π),pε[0,π),a,b,c are the radii, d is the center, and R isthe 3D rotation matrix of the ellipsoid.

According to an embodiment of the present disclosure, anellipse-segmenting method utilizes geodesic energy. If the contourC(t)=ε(t) is an ellipse with parameters λ_(i), the energy of the contouris:E(ε)=∫₀ ¹Φ∥ε_(p) ∥dp,where Φ is a weighting function, which is usually designed to slow downthe shrinkage of the contour at high image gradients.

Taking the derivative of the energy with respect to an independent timevariable t yields: $\begin{matrix}\begin{matrix}{{\int_{0}^{1}{{\frac{\mathbb{d}}{\mathbb{d}t}\left\lbrack {\Phi{\varepsilon_{p}}} \right\rbrack}\quad{\mathbb{d}p}}} = {{\int_{0}^{1}{\Phi_{t}{\varepsilon_{p}}\quad{\mathbb{d}p}}} + {\int_{0}^{1}{\Phi\frac{\mathbb{d}}{\mathbb{d}t}\sqrt{\left\langle {\varepsilon_{p},\varepsilon_{p}} \right\rangle}\quad{\mathbb{d}p}}}}} \\{= {{\int_{0}^{1}{\Phi_{t}{\varepsilon_{p}}\quad{\mathbb{d}p}}} + {\int_{0}^{1}{\Phi\left\langle {\varepsilon_{p\quad t},\underset{\underset{T}{︸}}{\varepsilon_{P}/{\varepsilon_{P}}}} \right\rangle\quad{\mathbb{d}p}}}}} \\{= {{\int_{0}^{1}{\Phi_{t}{\varepsilon_{p}}\quad{\mathbb{d}p}}} - {\int_{0}^{1}{\Phi\left\langle {\varepsilon_{t},T_{p}} \right\rangle\quad{\mathbb{d}p}}}}}\end{matrix} & \quad \\{\frac{\partial\lambda_{i}}{\partial t} = {{\int_{0}^{1}{\left\langle {{\nabla\Phi},\frac{\partial\varepsilon}{\partial\lambda^{i}}} \right\rangle{\varepsilon_{p}}\quad{\mathbb{d}p}}} - {\int_{0}^{1}{\Phi\left\langle {\frac{\partial\varepsilon}{\partial\lambda^{i}},T_{p}} \right\rangle{\mathbb{d}p}}}}} & (12)\end{matrix}$

-   -   where the T(p) is the ellipse tangent vector and T_(p) its        derivative (which is then in the normal direction of the        ellipse) along the ellipse, and the last line is the geodesic        evolution equation for the ellipse parameters λ_(i).

Taking a closer look at this flow, one can note the similarity to thegeodesic contour flow, where the second term is the curvature termweighted by the conformal factor Φ, and the first term is the gradientof the conformal factor that pulls the contour back to the realboundary. In the above equation though, the integration around theellipse provides a significant increase in robustness of the flow,allowing the contour to escape from local minima more easily, incontrast to a generic contour with geodesic energy.

A separate ellipse is utilized for geodesic joint segmentation andregistration. It is assumed that a common contour/ellipse E lives on anindependent domain, and let ε¹=g₁(ε), ε²=g₂(ε) (for simplicity inillustration, m=2 ellipses). $\begin{matrix}\begin{matrix}{{{E\left( {\varepsilon,g_{1},g_{2}} \right)} = {{\int_{\varepsilon^{1}}{{\Phi^{1}\left( x_{1} \right)}{{\varepsilon^{1}}_{p}}{\mathbb{d}p}}} + {\int_{\varepsilon^{2}}{{\Phi^{2}\left( x_{2} \right)}{\varepsilon_{p}^{2}}{\mathbb{d}p}}}}}\quad} \\{= {{\int_{0}^{1}{{\Phi^{1}\left( {g_{1}(x)} \right)}{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}{\mathbb{d}p}}} +}} \\{\int_{0}^{1}{{\Phi^{2}\left( {g_{2}(x)} \right)}{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}{\mathbb{d}p}}}\end{matrix} & (13)\end{matrix}$The evolution of the ellipse E is given by: $\begin{matrix}{\frac{\partial E}{\partial t} = {{\int_{0}^{1}{{\frac{\mathbb{d}}{\mathbb{d}t}\left\lbrack {{\Phi^{1}\left( {g_{1}(x)} \right)}{\varepsilon_{p}^{1}}} \right\rbrack}{\mathbb{d}p}}} + {\int_{0}^{1}{{\frac{\mathbb{d}}{\mathbb{d}t}\left\lbrack {{\Phi^{2}\left( {g_{2}(x)} \right)}{\varepsilon_{p}^{2}}} \right\rbrack}{\mathbb{d}p}}}}} \\{= {{\int_{0}^{1}{{\Phi_{t}^{1}\left( {g_{1}(x)} \right)}{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}{\mathbb{d}p}}} + {\int_{0}^{1}{{\Phi_{t}^{2}\left( {g_{2}(x)} \right)}{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}{\mathbb{d}p}}} +}} \\{{\left. {{\Phi^{1}\left( {g_{1}(x)} \right)}\frac{\mathbb{d}}{\mathbb{d}t}\sqrt{\left\langle {{g_{1}^{\prime}\left( \varepsilon_{p} \right)},{g_{1}^{\prime}\left( \varepsilon_{p} \right)}} \right\rangle}} \right){\mathbb{d}p}} +} \\{\left. {{\Phi^{2}\left( {g_{2}(x)} \right)}\frac{\mathbb{d}}{\mathbb{d}t}\sqrt{\left\langle {{g_{2}^{\prime}\left( \varepsilon_{p} \right)},{g_{2}^{\prime}\left( \varepsilon_{p} \right)}} \right\rangle}} \right){\mathbb{d}p}} \\{= {{\int_{0}^{1}{{\Phi_{t}^{1}\left( {g_{1}(x)} \right)}{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}}} + {{\Phi_{t}^{2}\left( {g_{2}(x)} \right)}{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}{\mathbb{d}p}} -}} \\{{\int_{0}^{1}{{\Phi^{1}\left( {g_{1}(x)} \right)}\left\langle {\left( {g_{1}^{\prime}\left( \varepsilon_{p} \right)} \right)_{t},\underset{\underset{T^{1}}{︸}}{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}/{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}}} \right\rangle{\mathbb{d}p}}} -} \\{\int_{0}^{1}{{\Phi^{2}\left( {g_{2}(x)} \right)}\left\langle {\left( {g_{2}^{\prime}\left( \varepsilon_{p} \right)} \right)_{t},\underset{\underset{T^{2}}{︸}}{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}/{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}}} \right\rangle{\mathbb{d}p}}} \\{= {{\int_{0}^{1}{{\Phi_{t}^{1}\left( {g_{1}(x)} \right)}{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}}} - {{\Phi^{1}\left( {g_{1}(x)} \right)}\left\langle {\left( {g_{1}^{\prime}(\varepsilon)} \right)_{t},T_{p}^{1}} \right\rangle{\mathbb{d}p}} +}} \\{{\int_{0}^{1}{{\Phi_{t}^{2}\left( {g_{2}(x)} \right)}{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}}} - {{\Phi^{2}\left( {g_{2}(x)} \right)}\left\langle {\left( {g_{2}^{\prime}(\varepsilon)} \right)_{t},T_{p}^{2}} \right\rangle{\mathbb{d}p}}}\end{matrix}$

Therefore, the evolution of the parameters of the ellipse λ_(i) aregiven by: $\begin{matrix}\begin{matrix}{\frac{\partial\lambda_{i}}{\partial t} = {{\int_{0}^{1}{\left\langle {{\nabla{\Phi^{1}\left( {g_{1}(x)} \right)}},\frac{\partial{g_{1}(\varepsilon)}}{\partial\lambda^{i}}} \right\rangle{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}}}\quad -}} \\{{\int_{0}^{1}{{\Phi^{1}\left( {g_{1}(x)} \right)}\left\langle {\frac{\partial{g_{1}(\varepsilon)}}{\partial\lambda^{i}},T_{p}^{1}} \right\rangle{\mathbb{d}p}}} +} \\{{\int_{0}^{1}{\left\langle {{\nabla{\Phi^{2}\left( {g_{2}(x)} \right)}},\frac{\partial{g_{2}(\varepsilon)}}{\partial\lambda^{i}}} \right\rangle{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}}}\quad -} \\{\int_{0}^{1}{{\Phi^{2}\left( {g_{2}(x)} \right)}\left\langle {\frac{\partial{g_{2}(\varepsilon)}}{\partial\lambda^{i}},T_{p}^{2}} \right\rangle{\mathbb{d}p}}}\end{matrix} & (14)\end{matrix}$

The curvature-like term, i.e. the term multiplied by the conformalfactor Φ^(i), is not of much use for an ellipse, since it is alwaysconvex. Therefore, such terms are eliminated. The influence of thegradient of Φ^(i) is extended by utilizing a gradient vector flow methodto prevent the flow in Eq. (14) from staying too local.

A geodesic ellipse registration evolution is employed. The energy forthe first registration g, is given by:E(g ₁)=∫₀ ¹Φ¹(x ₁)∥ε_(p) ¹ ∥dp=∫ ₀ ¹Φ¹(g ₁(x))∥g ₁′(ε_(p))∥dp  (15)

Then the evolution of the first registration g₁ can be obtained asfollows: $\begin{matrix}{{\frac{\partial\left( g_{1} \right)_{i}}{\partial t} = {\frac{\partial E}{\partial\left( g_{1} \right)_{i}} = {\int_{0}^{1}{\left\langle {{\nabla{\Phi^{1}\left( {g_{1}(x)} \right)}},\frac{\partial{g_{1}(\varepsilon)}}{\partial w_{i}}} \right\rangle{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}\quad{\mathbb{d}p}}}}},} & (16)\end{matrix}$and similarly the evolution of the second registration g₂$\begin{matrix}{\frac{\partial\left( g_{2} \right)_{i}}{\partial t} = {\frac{\partial E}{\partial\left( g_{2} \right)_{i}} = {\int_{0}^{1}{\left\langle {{\nabla{\Phi^{2}\left( {g_{2}(x)} \right)}},\frac{\partial{g_{2}(\varepsilon)}}{\partial w_{i}}} \right\rangle{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}\quad{{\mathbb{d}p}.}}}}} & (17)\end{matrix}$

The combined region and boundary based ellipse flow utilizes both theregion and boundary information in segmentation and registration ofstructures in medical images. These images have a high variability inthe condition of the border and intensity information of the targetstructure. A weighted combination of the region and geodesic flows areused to obtain the resulting update equations for both the registrationand the segmentation of the ellipse as follows: $\begin{matrix}\begin{matrix}{\frac{\partial\left( g_{k} \right)_{i}}{\partial t} = \frac{\partial E}{\partial\left( g_{k} \right)_{i}}} \\{= {\int_{0}^{1}\left\langle {\left\lbrack {{\alpha{\nabla{\Phi^{k}\left( {g_{k}(x)} \right)}}} + {\left( {1 - \alpha} \right){f^{k}\left( {g_{k}(x)} \right)}N}} \right\rbrack,\frac{\partial{g_{k}(x)}}{\partial w_{i}}} \right\rangle}} \\{{{{g_{k}^{\prime}\left( \varepsilon_{p} \right)}}\quad{\mathbb{d}p}},}\end{matrix} & (18) \\\begin{matrix}{\frac{\partial\lambda_{i}}{\partial t} = {\int_{0}^{1}\left\langle {\frac{\partial{g_{1}(\varepsilon)}}{\partial\lambda_{i}},\left\lbrack {{\alpha{\nabla{\Phi^{1}\left( {g_{1}(x)} \right)}}} + {\left( {1 - \alpha} \right){f^{1}\left( {g_{1}(x)} \right)}N}} \right\rbrack} \right\rangle}} \\{{{{g_{1}^{\prime}\left( \varepsilon_{p} \right)}}\quad{\mathbb{d}p}} +} \\{\int_{0}^{1}\left\langle {\frac{\partial{g_{2}(\varepsilon)}}{\partial\lambda_{i}},\left\lbrack {{\alpha{\nabla{\Phi^{2}\left( {g_{2}(x)} \right)}}} + {\left( {1 - \alpha} \right){f^{2}\left( {g_{2}(x)} \right)}N}} \right\rbrack} \right\rangle} \\{{{g_{2}^{\prime}\left( \varepsilon_{p} \right)}}\quad{\mathbb{d}p}}\end{matrix} & (19)\end{matrix}$for the ith parameter of both the registration and segmentation i=1, . .. , 5, and for the kth transformation to image k, (k=1,2 here forsimplicity). The flows in Eqs. (18-19) were tested on a set of threesynthetic objects, each subject to different rotation, translation andscaling. FIGS. 4A, 4D and 4G show an initialization of ellipses, e.g.,402 and 403. FIGS. 4B, 4E and 4H shown a respective interim step. FIGS.4C, 4F and 4I show a respective converged result. Each of FIGS. 4A-4Idepicts a synthetic object to be detected, e.g., 401, a seed position402 and an ellipse, e.g., 403. An equal weighting from both geodesic andregion-based terms is used, where α=0.5. Note that these examples alsodemonstrates that an initialization that is only good for some imagesstill leads to a stable result on all images.

General Contour with Non-Rigid Registration

The non-rigid deformation corresponds to a registration in infinitedimensions without using shape priors. If required by the applicationthough, the shape priors can be incorporated into the framework. Thejoint segmentation and registration accounts for general problems ofregistration among anatomical structures defined by a deformation fieldor a displacement field between target regions.

Referring to FIG. 5A, a joint non-rigid registration and segmentationmethod includes initializing a seed on a slice of an image volume 501,for example, manually selecting the seed. Based on the initialized seed,a level set function is initialized 502 that is a signed distancefunction, negative inside and positive outside the surface, which isconverted into a 3D surface seed by extending to two neighboring slices.A non-rigid registration is initialized 503, which is a deformationvector field to a zero vector field over the image domain. The surface'slevel set function and deformation vector field are updated 504, and theupdated surface, either its signed distance function over the domain orthe surface that is the zero-level set of the signed distance function,and the updated deformation field U(TotalTime), which maps the surfacein one of the images onto another image, are output 505. FIG. 6illustrates a point 601 in a first domain 602 being transformed ormapped into a second domain 603.

Given two image volumes, the transformation that deforms one of theimages to the other one may be denoted by: {circumflex over(x)}=T(x)=x+u(x), where u(x) is the displacement vector field defined asu:Ω→

^(n) (n=2 or 3).

The general region-based energy functional for a joint segmentation andregistration purpose can be rewritten as follows:E(S,u)=∫_(S) _(in) f _(in)(x)dx+∫ _(S) _(out) f _(out)(x)dx+∫ _(S) dA∫_(Ŝ) _(in) {circumflex over (f)} _(in)({circumflex over(x)})d{circumflex over (x)}+∫ _(Ŝ) _(out) {circumflex over (f)}_(out)({circumflex over (x)})d{circumflex over (x)}+∫ _(Ω) ∥Du∥ ²dx  (20)

-   -   where the last two terms represent the regularization on the        unknown surface S and vector field u,{circumflex over        (x)}=x+u(x), and f_(in), f_(out), {circumflex over (f)}_(in) and        {circumflex over (f)}_(out) are as before.

The solutions to the minimization problems may be given by:$\begin{matrix}{{\overset{\sim}{S} = {\arg\quad{\min\limits_{S}{E\left( {S,u} \right)}}}},{{{and}\quad\overset{\sim}{u}} = {\arg\quad{\min\limits_{u}{E\left( {S,u} \right)}}}}} & (21)\end{matrix}$

Let us reformulate the problem in terms of a level set function Φ:Ω→

that represents S as its zero level set. Let H denote a Heavisidefunction, then E in Eq. (20) can be re-written as:E(Φ,u)=∫_(Ω) f _(in)(x)H(Φ(x))dx+∫ _(Ω) f _(out)(x)(1−H(Φ(x)))dx+∫_({circumflex over (Ω)}) {circumflex over (f)} _(in)({circumflex over(x)})H({circumflex over (Φ)}({circumflex over (x)}))d{circumflex over(x)}+∫ _({circumflex over (Ω)}) {circumflex over (f)} _(out)({circumflexover (x)})(1−H({circumflex over (Φ)}({circumflex over(x)})))d{circumflex over (x)}+∫ _(Ω)δ(Φ)|∇Φ|dx+∫ _(Ω) ∥Du(x)∥² dx  (21)

-   -   where {circumflex over (Φ)}(x+u)=Φ(x) and δ(x)=dH(x)/dx in the        sense of distributions. Regularized versions of H(x) and δ(x)        may be used, for example, according a Heaviside function and        delta function, respectively. The surface evolution in this case        is given by: $\begin{matrix}        {\frac{\partial\Phi}{\partial t} = {{{\delta(\Phi)}{f(x)}} + {{\delta\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}{\hat{f}\left( {x + {u(x)}} \right)}} + {\kappa.}}} & (22)        \end{matrix}$

For registration evolution, the only part of the energy functional inEq. (21) that is taken into account is: $\begin{matrix}{{E(u)} = {{\int_{\Omega}{\underset{\underset{F_{i\quad n}{({x + {u{(x)}}})}}{︸}}{{{\hat{f}}_{i\quad n}\left( {x + {u(x)}} \right)}{H\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}}\quad{\mathbb{d}x}}} + {\int_{\Omega}{\underset{\underset{F_{out}{({x + {u{(x)}}})}}{︸}}{{{\hat{f}}_{out}\left( {x + {u(x)}} \right)}\left( {1 - {H\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}} \right)}\quad{\mathbb{d}x}}} + {\alpha{\int_{\Omega}{{{Du}}^{2}\quad{\mathbb{d}x}}}}}} & (23)\end{matrix}$

The gradient of E with respect to u is${\frac{\partial E}{\partial u} = {{{- {\nabla u}}\quad{F_{in}(u)}} - {{\nabla u}\quad{F_{out}(u)}} + {\alpha\quad\Delta\quad u}}},$and we solve it by a gradient descent PDE (over domain Ω) whose steadystate solution gives the minimizer displacement field u that varies overspace: $\begin{matrix}{\frac{\partial{u\left( {x,t} \right)}}{\partial t} = {{{\nabla u}{{\hat{f}}_{i\quad n}\left( {x + u} \right)}{H\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}} + {{\nabla u}{{\hat{f}}_{out}\left( {x + u} \right)}\left( {1 - {H\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}} \right)} + {{\hat{f}\left( {x + u} \right)}{\delta\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}{\nabla{\hat{\Phi}\left( {x + u} \right)}}} + {{\alpha\Delta}\quad{u(x)}}}} & \quad \\{{u\left( {x,0} \right)} = {u_{o}(x)}} & (24)\end{matrix}$

Due to a narrowband level set implementation, the method is effectivelysolving for u(x) on a band around the surface.

Referring to FIG. 5B, the following method may be used for updating thesurface's level set function and the deformation vector field. ΔU(X,t) =0, for all X in the image domain, ΔΦ(t) = 0, for all X in the imagedomain where X is a coordinate, e.g., (x,y) for 2D or (x,y,z) for 3D ForX in the image volume domain (506) //Accumulate gradients for surfaceupdate (508) ΔΦ(X,t)+= right portion of Eq.22 //Accumulate gradients fordeformation field update (508) ΔU(X,t)+= right portion of Eq.24 End ForX in the image volume domain //Update the segmentation of the surface(509) Φ(X,t + Δt) = Φ(X,t) + ΔtΔΦ(X,t) //Update the deformation field(510) U(X,t + Δt) = U(X,t) + ΔtΔU(X,t) //Discretized Eq.24 End

If as a special case, a piecewise (p-w) constant model known as“Chan-Vese” type flows is chosen for the target regions that are to besegmented and registered in images I and Î, the region based term in theenergy, and the gradient terms can be given by: $\begin{matrix}{\hat{f} = {{\hat{f}}_{i\quad n} - {\hat{f}}_{out}}} \\{{= {2\left( {{\hat{m}}_{i\quad n} - {\hat{m}}_{out}} \right)\left( {\frac{{\hat{m}}_{i\quad n} + {\hat{m}}_{out}}{2} - {\hat{I}\left( {x + {u(x)}} \right)}} \right)}},}\end{matrix}$ ∇uf̂_(i  n) = 2(Î(x + u) − m̂_(i  n))∇Î(x + u(x))∇uf̂_(out) = 2(Î(x + u) − m̂_(out))∇Î(x + u(x)).

These expressions can be inserted into Eq. (24) to obtain the PDE forevolution of non-rigid registration field for the p-w constant regionmodel: $\begin{matrix}{\frac{\partial{u\left( {x,t} \right)}}{\partial t} = {{{- {\left( {{\hat{m}}_{i\quad n} - {\hat{m}}_{out}} \right)\left\lbrack {\frac{\left( {{\hat{m}}_{i\quad n} + {\hat{m}}_{out}} \right)}{2} - {\hat{I}\left( {x + u} \right)}} \right\rbrack}}\delta\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right){\nabla{\hat{\Phi}\left( {x + u} \right)}}} - {\left( {{\hat{I}\left( {x + u} \right)} - {\hat{m}}_{i\quad n}} \right){\nabla{\hat{I}\left( {x + {u(x)}} \right)}}{H\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}} - {\left( {{\hat{I}\left( {x + u} \right)} - {\hat{m}}_{out}} \right){\nabla{\hat{I}\left( {x + {u(x)}} \right)}}\left( {1 - {H\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right)}} \right)} + {\frac{\alpha}{2}\Delta\quad{u(x)}}}} & \quad \\{{{u\left( {x,0} \right)} = {{u_{o}(x)} = 0}},} & (25)\end{matrix}$

-   -   where a zero vector field initialization is adequate in solving        the PDE without any prior knowledge of the true vector field u.

In some segmentation applications, the basic approach of thresholdinghas proven to be useful. Instead of using means inside and outside thesurface, one can convert such a basic “Chan-Vese” flow to a“thresholding” flow. This is also equivalent to region growing toseparate the intensity inside the growing surface from the outside bythe given threshold. For this purpose the following speed function maybe used:{circumflex over (f)}=({circumflex over (m)} _(in) −{circumflex over(m)} _(out))(T−Î(x+u(x)))  (26)

-   -   where $\frac{{\hat{m}}_{i\quad n} + {\hat{m}}_{out}}{2}$        quantity in “Chan-Vese” flow is replaced by an arbitrary        threshold T. For this speed function the following PDE may be        used for updating the vector field: $\begin{matrix}        \begin{matrix}        {\frac{\partial{u\left( {x,t} \right)}}{\partial t} = {{- \left( {{\hat{m}}_{i\quad n} - {\hat{m}}_{out}} \right)}\left( {T - {\hat{I}\left( {x + u} \right)}} \right)}} \\        {{\delta\left( {\hat{\Phi}\left( {x + {u(x)}} \right)} \right){\nabla{\hat{\Phi}\left( {x + u} \right)}}} + {\frac{\alpha}{2}\Delta\quad{u(x)}}}        \end{matrix} & \quad \\        {{{u\left( {x,0} \right)} = {{u_{o}(x)} = 0}},} & (27)        \end{matrix}$        with the boundary term only.

Having described embodiments for a system and method for a jointsegmentation and registration of images for object detection, it isnoted that modifications and variations can be made by persons skilledin the art in light of the above teachings. It is therefore to beunderstood that changes may be made in the particular embodiments of theinvention disclosed which are within the scope and spirit of theinvention as defined by the appended claims. Having thus described theinvention with the details and particularity required by the patentlaws, what is claimed and desired protected by Letters Patent is setforth in the appended claims.

1. A method for joint segmentation and registration of a contour fordetecting an object in image data comprising: initializing asegmentation; initializing a registration of the segmentation; updatingthe segmentation and the registration; and outputting an updatedsegmentation and an updated registration, wherein the updatedsegmentation and updated registration delineates the object in the imagedata.
 2. The method of claim 1, wherein the segmentation is an ellipseand the registration is rigid.
 3. The method of claim 2, whereinupdating the segmentation and registration further comprises:accumulating gradients for the segmentation evolution for a plurality ofimages in the image data and a plurality of samples on the ellipse;updating the registration for each of the plurality of images; andupdating the segmentation given the updated registration.
 4. The methodof claim 2, further comprising determining parameters of each rigidregistration including a translation vector, a rotation matrix and anon-uniform scale matrix.
 5. The method of claim 1, wherein thesegmentation is a surface of a three-dimensional structure and theregistration is non-rigid.
 6. The method of claim 5, wherein thenon-rigid registration is a deformation vector field over the surface ofthe three-dimensional object.
 7. The method of claim 6, wherein updatingthe segmentation and registration further comprises: accumulatinggradients for surface evolution for each coordinate in the image volume;and accumulating gradients for the non-rigid registration evolution foreach coordinate in the image volume.
 8. The method claim 6, wherein thedeformation vector field maps the surface in a first image onto thesurface of in a second image.
 9. The method of claim 5, wherein thethree-dimensional surface is determined from an initialized surface seedof the segmentation in a first image slice of the image data and twoneighboring slices in the image data in which the surface seed isextended.
 10. A program storage device readable by machine, tangiblyembodying a program of instructions executable by the machine to performmethod steps for joint segmentation and registration of a contour fordetecting an object in image data, the method steps comprising:initializing a segmentation; initializing a registration of thesegmentation; updating the segmentation and the registration; andoutputting an updated segmentation and an updated registration, whereinthe updated segmentation and updated registration delineates the objectin the image data.
 11. The method of claim 10, wherein the segmentationis an ellipse and the registration is rigid.
 12. The method of claim 11,wherein updating the segmentation and registration further comprises:accumulating gradients for the segmentation evolution for a plurality ofimages in the image data and a plurality of samples on the ellipse;updating the registration for each of the plurality of images; andupdating the segmentation given the updated registration.
 13. The methodof claim 10, wherein the segmentation is a surface of athree-dimensional structure and the registration is non-rigid.
 14. Themethod of claim 13, wherein the non-rigid registration is a deformationvector field over the surface of the three-dimensional object.
 15. Themethod of claim 14, wherein updating the segmentation and registrationfurther comprises: accumulating gradients for surface evolution for eachcoordinate in the image volume; and accumulating gradients for thenon-rigid registration evolution for each coordinate in the imagevolume.
 16. A method for joint segmentation and registration of acontour for detecting an object in an image volume comprising:initializing a seed in a first slice of the image volume; initializing asigned distance function; converting the distance function into athree-dimensional surface seed; initializing a non-rigid registrationfor the surface; updating the distance function and the non-rigidregistration; and outputting an updated surface and distance functionmapping the surface in the first slice to the surface in a second slice,the surface delineating the object in the image volume.
 17. The methodof claim 16, wherein updating the segmentation and registration furthercomprises: accumulating gradients for surface evolution for eachcoordinate in the image volume; and accumulating gradients for thenon-rigid registration evolution for each coordinate in the imagevolume.
 18. The method of claim 16, wherein the surface is determinedfrom an initialized surface seed of the segmentation in the first imageslice of the image data and two neighboring slices in the image data inwhich the surface seed is extended.